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The role of frustration and interaction strength on the half-filled Hubbard model is studied on 
the square lattice with nearest and next-nearest neighbour hoppings t and t' using the Variational 
Cluster Approximation (VGA). At half- filling, we find two phases with long-range antiferromagnetic 
(AF) order: the usual Neel phase, stable at small frustration t' /t, and the so-called coUinear (or 
super-antiferromagnet) phase with ordering wave- vector (tt, 0) or (0, tt), stable for large frustration. 
These are separated by a phase with no detectable long-range magnetic order. We also find the 
d-wave superconducting (SC) phase {d^2_y2), which is favoured by frustration if it is not too large. 
Intriguingly, there is a broad region of coexistence where both AF and SC order parameters have non- 
zero values. In addition, the physics of the metal-insulator transition in the normal state is analyzed. 
The results obtained with the help of the VGA method are compared with the large- C/ expansion 
of the Hubbard model and known results for the frustrated J1-J2 Heisenberg model. These results 
are relevant for pressure studies of undoped parents of the high-temperature superconductors: we 
predict that an insulator to d-wave SG transition may appear under pressure. 

PACS numbers: 71.10.Fd, 74.72.-h, 71.30.+h, 75.10.Jm 



I. INTRODUCTION 

The subject of frustration^ in quantum magnetic sys- 
tems has received increased attention in recent years, fu- 
elled in part by the discovery^ of high-temperature su- 
perconductivity in the doped cuprates. It is beheved to 
play a key role in a number of recently observed phenom- 
ena, such as the large anomalous Hall effect in ferromag- 
netic pyrochlores^, the unconventional superconductivity 
in water substituted sodium cobaltate Naa;Co02, which 
is composed of triangular sheets of Co atoms'*, the inter- 
play between magnetism and unconventional supercon- 
ductivity in organic layered compounds of the k-BEDT 
family^'^ or the interaction between electric and magnetic 
properties in multiferroic materials^. 

The issue of frustration has been studied mostly on two 
classes of theoretical models: spin Hamiltonians, such 
as the J1-J2 Heisenberg model discussed below, and toy 
dimer models, the latter inspired by P.W. Anderson's 
proposal* of the resonating valence-bond (RVB) state as 
a possible explanation for high-Tc superconductivity. 

One can view spin Hamiltonians as the large interac- 
tion, U, limit of the Hubbard model. It is thus of interest 
to study the effect of both interaction and frustration on 
the phase diagram. In this work, we study systemati- 
cally the frustrated Hubbard model at half-filling on a 
square lattice with nearest t and next-nearest neighbour 
t' hoppings, described by the Hamiltonian: 

^ = i E ^^J- + II ^^^^^ "'T"U, (1) 

<i,j> <€.i,j:^ i 

where c|, Ci are the electron creation and annihilation 
operators and riic is the particle number operator on 
site i. The interaction is represented by U . Next-nearest 



neighbor hopping t' introduces frustration since, from 
a weak coupling point of view, it produces deviations 
from perfect nesting, and, from a strong-coupling point 
of view, it leads to an effective antiferromagnetic su- 
perexchange interaction J2 that opposes the tendency of 
next-nearest neighbors to order ferromagnetically when 
nearest-neighbor superexchange Ji is antiferromagnetic. 

Our study of the phase diagram as a function of 
U /t and t' /t can also be understood as a study of the 
generalized zero-temperature phase diagram for high- 
temperature superconductors illustrated in Fig. 1. The 
thin parallelepiped represents schematically the region of 
parameter space where families of high-temperature su- 
perconductors appear. We are studying the zero-doping 
plane (5 = 0, where one normally encounters the insu- 
lating antiferromagnetic parents of high-temperature su- 
perconductors. We will see that d-wave superconductiv- 
ity can also occur in this plane, so high-pressure studies 
might conceivably lead to the observation of d-wave su- 
perconductivity even at half-filling, provided the on-site 
interaction U /t is not too large. The generalized phase 
diagram also leads to insights into the nature of d-wave 
superconductivity, as we will see. 

We use a quantum cluster approach, the so-called vari- 
ational cluster approximation (VCA, sometimes referred 
to as the 'variational cluster perturbation theory')^. This 
method has already been used successfully for the high- 
temperature superconductors^''^^^. Other quantum clus- 
ter methods that are extensions of dynamical mean-field 
theory (DMFT)^^, such as cellular dynamical mean-field 
theory and dynamical cluster approximation^^, have 
yielded comparable results^^'^^ for the same problem. 
On the anisotropic triangular lattice at half-filling, both 
VCA** and CDMFT*^ give a phase diagram that is in re- 
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<i, j> and respectively) : 



H 



Ji ^ Si ■ Sj + J2 ^ Si • Sk 



(2) 



FIG. 1: (Color online) Schematic generalized zero- 
temperature parameter space for the high-temperature super- 
conductors. Horizontal axis 5 represents doping, vertical axis 
U/t interaction strength, and third direction t' /t frustration. 
The parallelepiped indicates the region of parameter space 
relevant for high-temperature superconductors. In this study, 
we consider the zero doping 5 = plane. Experimentally, U/t 
can be varied by pressure. The inset shows the definitions of 
t and t' on the square lattice. 



markable agreement with that of layered BEDT organic 
superconductors. 

Wherever possible, we will make connections with 
earlier numerical studies on the half-filled square- 
lattice Hubbard model^"'^^, and with work on spin- 
Hamiltonians. 

This paper is organized as follows. First, the role 
of frustration in quantum magnets is illustrated in Sec- 
tion II on the well-studied example of the J1-J2 Heisen- 
berg spin model. The connection to the Hubbard model 
is then made by virtue of the large-[/ expansion in Sec- 
tion HI. The framework of the variational cluster ap- 
proximation (VGA) used in this work is briefly described 
in Section IV. The resulting magnetic phase diagram of 
the Hubbard model is then presented in Section V, with a 
separate subsection devoted to the analysis of the metal- 
insulator transition. The main result of this work, where 
d-wave superconductivity and magnetism compete and 
even coexist, is described in Section VI. We conclude by 
discussing the obtained phase diagram of the frustrated 
Hubbard model and draw comparison with other known 
results in Section VII. 



II. A REFERENCE POINT: J1-J2 HEISENBERG 
MODEL ON THE SQUARE LATTICE 

One of the earliest studied models that exhibits frus- 
tration is the so-called J1-J2 Heisenberg model, which 
contains antiferromagnetic (AF) spin-spin interaction be- 
tween nearest and next-nearest neighbours (denoted by 



<4J> 



Albeit simple in appearance, this model captures a num- 
ber of important features common to a large class of frus- 
trated quantum magnets. 

Classically, the ground state of the model can be de- 
rived by considering the Fourier transform of the spin 
coupling J{q), which on the square lattice with next- 
nearest neighbour spin interaction takes the following 
form: 

J{q) ~ 2Ji{cosqj; + cosQy) + 4J2 cosq^ cosqy. (3) 

The classical ground state should minimize this cou- 
pling, leading to two possible solutions: the Neel state 
(referred to as AFl in the following) with the order- 
ing wave-vector Q = (tt, tt) for the range of parameters 
J2/J1 < 0.5, and the so-called super-antiferromagnetic 
phase with Q = (tt, 0) or (0, tt) (referred to as AF2), 
realised for J2/J1 > 0.5. Although non-coUinear spin 
states with the same classical ground state energy can 
also be realized, it has been shown^^ that thermal or 
quantum fluctuations will favour the states that have 
coUinear magnetization. 

The effect of quantum fluctuations becomes especially 
important around the quantum critical point J2/ Ji — 0.5 
where the classical ground state is highly degenerate. 
The large-S' analysis shows^^ that even to the lowest 
order in l/S, zero-temperature quantum corrections to 
the sublattice magnetization diverge at the critical point, 
pointing to the existence of a quantum disordered phase. 
The nature of such a phase can be captured by dimer cov- 
ering of the lattice, which is a caricature for the singlet 
pairings (i.e. valence bonds) of nearest-neighbour spins. 

A wide literature^'' exists on the subject of spin ro- 
tationally invariant dimer order in frustrated quantum 
magnets. It is generally believed^^ that in the case of a 
square lattice, the dimer phase exhibits long-range order 
in the dimer-dimer correlation functions, leading to the 
notion of the 'valence bond solid' (VBS), as opposed to 
the original RVB phase of Anderson^ which is supposed 
to have only short-range order and gapped collective ex- 
citations. We shall touch upon this subject in Section V, 
although this study will be primarily concerned with the 
magnetic broken-symmetry phases. 



III. LARGE-[7 EXPANSION OF THE t - t' - U 
HUBBARD MODEL 

In order to get insight into the physics of the frustrated 
Hubbard model, we shall first consider its large-?/ expan- 
sion. Whereas the procedure for obtaining the low-energy 
Heisenberg Hamiltonian from the conventional Hubbard 
model with nearest-neighbour interaction is a textbook 
example^^, the presence of next-nearest neighbour terms 
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FIG. 2: The ring exchange contributions to the large- [/ ex- 
pansion of the t — t' — U Hubbard model from Ref. 26. The 
empty circle in case (c) denotes that the central site does not 
participate in the plaquette spin-exchange term. 

and calculation to order leads non-trivial next- 

nearest-neighbor and ring-exchange terms. We exploit 
here the coefficients of the expansion that have been ob- 
tained by Delannoy et aZ.^^'^"^ by means of the canoni- 
cal transformation approach. The resulting effective spin 
Hamiltonian can be written as follows: 

H — Ji ^ Si ■ Sj + J2 ^ Si • Sfc 

+ {rurg exchange terms} (4) 
where the coefficients J; arc given by 

The relevant ring exchange terms^^ arc defined on the 
plaquettes depicted in Fig. 2 with the corresponding an- 
alytical expressions given by: 

= 20Jtt'Pi{i,j,k,l) - Ju'P2{t,j,k,l) 
i7(b) = 20 Jw Pi{i, J, k, I) + Ju'P2ii,J,k,l) (6) 

-^(c) = 80■^Fl(^,J■ fc,0 - Ju'P2{t,j,k,l) 

where Jtt' = At^t'"^ /U^ and the following notations have 
been used following Ref. 26: 

Pi{i,j,kj) ^ (s, •s,)(Sfe.so + (s, •so(Sfc.s,) 

- (S..Sfe)(S,.Sz) (7) 
P2{i,j,k,l) = {Si ■ Sj + Si -Sk + Si -Si 

+ Sj-Sk + Sj-Si + Sk-Si} (8) 

Evaluating the classical ground state energies of the 
Hamiltonian Eqs. (4-7) that corresponds to the two pos- 
sible ordering wave- vectors Qi = (tt, tt) and Q2 = (ti", 0), 
yields the following result: 

2^ ( 12i'2 t'2\ ^ ^ 

i?(...)-i?K0) = ^(^-l + -^-^ + 2^j (9) 
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FIG. 3: (Color online) Classical phase diagram of the Hub- 
bard model with next-nearest neighbour hopping t' based on 
the large-C/ expansion analysis. The usual AF phase (AFl) 
with the ordering vector (tt, tt) is observed at small ratios of 
t'/t, whereas a so-called super-antiferromagnetic phase (0, tt) 
(AF2, shaded region) is realized at large frustration t' jt. Note 
that because of the nature of the large-(7 expansion, this phase 
diagram is expected to be inaccurate at low U values. 



The corresponding classical phase diagram that follows 
from this is shown in Fig. 3. Note that in the large of 
U limit, it follows from Eq. (9) that the criterion for the 
(tt, 0) phase to have lower ground state energy is given by 
t' jt > I/V2, i.e. J2/J1 > 0.5, which coincides with the 
classical criterion obtained earlier for the J1-J2 Heisen- 
berg model. 

Several comments on the phase diagram in Fig. 3 are 
due. First of all, the above analysis is based on the large- 
U expansion of the Hubbard model where electrons are 
localized. It should not be taken seriously for small U 
values where electrons can be delocalized even at half- 
filling. In particular, the AF2 (tt, 0) phase for a large 
range of t' /t below U < 2.8* is an artifact. Secondly, 
the role of thermal and quantum fluctuations has been 
completely neglected in the above classical analysis. The 
latter are going to lower the ground state energies of the 
two AF phases but, more importantly, they may favour 
a different type of order with no broken rotational spin 
symmetry, such as the valence bond solid or the RVB 
spin liquid state mentioned earlier in relation with the 
J1-J2 Heiscnberg model. 

Below, we take into account the effect of short-range 
quantum fluctuations at zero temperature as well as the 
possibility for electrons to delocalize. In order to achieve 
this we perform VGA calculations of the Hubbard model 
with nearest and next-nearest neighbour hopping and ob- 
tain the quantum analogue of the phase diagram shown 
in Fig. 3, including the possibility of d-wave supercon- 
ductivity. These results are presented in Section V. 
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IV. VARIATIONAL CLUSTER APPROACH 

Despite the apparent simpUcity of the Hubbard model, 
its phase diagram is extremely rich in physical phenom- 
ena, such as antiferromagnetism, incommensurate spin- 
density wave and d-wave superconductivity. Analytical 
progress is severely hindered by the fact that the model 
does not have a small parameter in the interesting regime, 
and, consequently, a number of numerical methods have 
been proposed to treat the Hubbard model. Among 
these, quantum cluster methods^* form a separate group, 
which seems to successfully capture many physical fea- 
tures of the model, including d-wave superconductivity 
that was studied^"'^^'^^'^^ in the context of the cuprates. 

In this work we have used the so-called variational clus- 
ter approximation (VGA), sometimes referred to as vari- 
ational cluster perturbation theory (VCPT) in the lit- 
erature. It is a special case of the self-energy functional 
approximation'^"'^^. The idea of this approach consists in 
expressing the grand canonical potential f2 of the model 
as a functional of the self-energy E: 



= Fm - Trln (-Gq ^ + S) 



(10) 



where Gq is the bare single-particle Green's function of 
the problem and F[T^] is the Legendre transform of the 
Luttinger-Ward functional, the latter being defined by 
the infinite sum of vacuum skeleton diagrams. 

The functional (10) can be proven to be stationary at 
the solution of the problem, i.e. where S is the physical 
self-energy of the Hubbard model. 



= 0. 



sol 



(11) 



The problem of finding the solution is then reduced to 
minimizing the functional f2[E] with respect to the self- 
energy S. Two problems stand in the way however. 
Firstly, the functional F^E] entering Eq. (10) is unknown 
and thus has to be approximated in some way. And sec- 
ondly, there is no easy practical way of varying the grand 
potential with respect to the self-energy. 

Potthoff suggested^° an elegant way around these 
problems by noting that since the functional F[T,] is a 
universal functional of the interaction part of the Hamil- 
tonian only (i.e. the last term in Eq. (1)), it can be 
obtained from the known (numerical) solution of a sim- 
pler reference system with the Hamiltonian H' defined 
on a partition of the infinite lattice into disjoint clusters, 
provided that the interaction term is kept the same as 
in the original Hamiltonian. For such a cluster partition, 
Eq. (10) can now be rewritten as 



- Trim 



(12) 



where the prime denotes the quantities defined on the 
cluster, to distinguish from those of the original problem. 
Gombining Eqs. (10) and (12), we finally obtain 



n[T,] = r^'p] +Trlni 



Trln(-Gn 



(13) 



Equation (13) is the central equation of the variational 
cluster approximation. The role of variational variables is 
played by some one-body parameters {h'} of the cluster 
Hamiltonian, so that one looks for a stationary solution 



an 

5/7 



SE dh' 



0. 



(14) 



It has been shown by Potthoff'^° that the VGA and an- 
other widely known method, cluster dynamical mean 
field theory (GDMFT), can both be formulated in the 
framework of the above self-energy formalism. The par- 
ticular advantage of the VGA is that it enables one to eas- 
ily study broken-symmetry phases for clusters of varying 
sizes. The Weiss fields h' are introduced into the clus- 
ter Hamiltonian and the potential Q is minimized with 
respect to it. It is important to stress that, in contrast 
with the usual mean-field theories, these Weiss fields are 
not mean fields, in a sense that the interaction part of 
the Hamiltonian is not factorized in any way and short- 
range correlations are treated exactly. The Weiss fields 
are introduced simply to allow for the possibility of a spe- 
cific long-range order, without ever imposing this order 
on the original Hamiltonian. 

In this work we have defined the cluster Hamiltonian 
with appropriate Weiss fields as follows: 



/C^tCx'i + h.c. 



-M J2 e''^"(-l)""x. - 51 "x,. + C/ 5] "xT^xi (15) 

x,cr X(T X 

Here as before, c^.^ is the electron creation operator at 
site X with spin cr, n^a is the particle number opera- 
tor, and M and Axx' are the Weiss fields corresponding 
to the antiferromagnetic order parameter (with ordering 
wave- vector Q) and to the superconducting order param- 
eters respectively. For singlet superconductivity we have 
Axx' — Ax'x- In particular, for d^2_y2 symmetry the 
Weiss field is defined as follows (e is a lattice vector): 



A, 



D, for e = ±2; 
—D, for e ~ ±y. 



(16) 



The corresponding order parameters, M and D, are given 
by the terms multiplying Af and D respectively in the 
Hamiltonian (15). 

In addition to the AF and SG Weiss fields, we also al- 
low the cluster chemical potential /x' to vary, to insure 
internal thermodynamic consistency'^^ of the calculation. 
Therefore, in all calculations reported in the present 
work, the cluster chemical potential fi' was treated as 
a variational parameter, along with symmetry-breaking 
Weiss fields, such as the staggered magnetization in the 
AF case. 

We solve the cluster problem using the Lanczos ex- 
act diagonalization technique, which enables one to find 
the ground state of the model at zero temperature. The 
cluster Green's function G|jj,(u;, k), defined for a pair of 
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FIG. 4: The tilings of the square lattice with the various clus- 
ters used in this work containing (a) 4 sites, (b) 6 sites, (c) 8 
sites. In this example the grey and white sites are inequiva- 
lent since the (7r,7r) AF order is possible. Note that in some 
cases, such as that of the 6-site cluster, a different tiling (d) 
must be chosen for a (0, vr) phase to become possible. 



cluster sites (a, b), was then evaluated using the so-called 
band Lanczos method^"^ that offers a significant compu- 
tational advantage over other approaches^*. The search 
for a stationary solution Eq. (14) was performed using a 
combination of the Newton-Raphson"^® and the conjugate 
gradient methods'^^. 

Since all the calculations are performed in the grand- 
canonical ensemble, the requirement on the filling (n) ~ 
1 is achieved by choosing the appropriate value of the 
lattice chemical potential ^. We note that this task is not 
trivial in the present study since the variation of both the 
cluster chemical potential /i' and the Weiss field (the SC 
D and AF M) tend to greatly influence the value of (n) . 
The appropriate value of the lattice chemical potential 
fi therefore had to be chosen at each point in the phase 
space of the parameters t, t' , U of the Hamiltonian (1) to 
guarantee that the system always remained at half-filling. 

In general the phase diagram will depend on the choice 
of the reference cluster system H' that is solved numer- 
ically to obtain the quantities entering the VGA equa- 
tion (13). The VGA solution becomes exact only in the 
thermodynamic limit of infinitely large cluster. In prac- 
tice, the typical cluster size is limited to a maximum of 
10-12 sites since the Hilbert space of the reference clus- 
ter Hamiltonian grows exponentially with the cluster size 
and so does the computational cost of the exact diagonal- 
ization algorithm. In this work we have studied clusters 
of 4, 6 and 8 sites, as depicted in Fig. 4. This appears 
sufficient to suggest what the result should look like in 
the thermodynamic limit. 



V. MAGNETISM AND MOTT PHYSICS IN 
THE FRUSTRATED HUBBARD MODEL 

A. Magnetic phase diagram 

When applying the VGA method to the frustrated 
Hubbard model, we have studied the possibilities for both 




FIG. 5: (Color online) The magnetic phase diagram of the 
t — t' — U Hubbard model ignoring the SC solution. There are 
two magnetic phases: (tt, tt), (tt, 0) and the paramagnetic re- 
gion (shaded area) where both order parameters vanish. The 
diagram was obtained with the VGA method using an 8-site 
cluster. The dashed line with filled data points shows the 
Hartree-Fock prediction from Ref. 37 for the transition be- 
tween the Neel (tt, tt) and the non-magnetic phases. The 
solid line (red online) with error-bars indicates the transition 
Uc{t') between the insulating (Mott-Iike) phase above and the 
metallic region below Uc when no magnetic order is allowed. 



long-range magnetic order and (d-wave) superconductiv- 
ity. In order to shed more light on the interplay between 
frustration and magnetism and to connect with existing 
studies on spin Hamiltonians summarized in Sec. II, we 
first report our results for purely magnetic phases, i.e. ig- 
noring the SG solution for the moment. The main results 
of this study can be summarized by the phase diagram in 
Fig. 5 where the horizontal axis is a measure of the frus- 
tration t' /t and the vertical axis the interaction strength 
U/t. 

We have looked for the same two AF phases that are 
predicted by both the J1-J2 Heisenberg model (Sec. II) 
and the large-?/ expansion of the Hubbard model (see 
Sec. HI), namely the usual Neel phase with the order- 
ing wave-vector Q = (tt, tt) and the so-called collinear 
order with Q ~ (tt, 0) (or equivalently, (0,7r)). The re- 
gions of stability of these two phases are shown on the 
phase diagram in Fig. 5 for the largest cluster studied 
(8 sites). The two magnetic phases are separated by 
a paramagnetic region (filled area in Fig. 5) where no 
non-zero value was found for either of the two order pa- 
rameters. We shall refer to this paramagnetic region as 
"disordered" although, strictly speaking, we cannot ex- 
clude the possibility of some other magnetic long-range 
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FIG. 6: Fermi surface for t' = 0.2t (left), t' = 0.6t (center) and 
t' — 0.8t (right). The change in topology (Lifshitz transition) 
occurs around t' = 0.71t. The best nesting vector Q is also 
shown. The Fermi surface for negative t' /t looks like the ones 
above if one translates the origin to (tt, tt). The phase diagram 
in other figures depends only on the absolute value of t' /t. 



order with, for example, incommensurate wavevector Q, 
which is not tractable by our method because of the finite 
cluster sizcs^^. 

As is clear from Fig. 5, for large U values the disordered 
region is centred around the critical value t'/t = 1/V2, 
confirming the predictions of the large-?/ expansion (c.f. 
Fig. 3). This region then becomes broader as U is low- 
ered, engulfing the whole of the phase diagram in the 
limit oi U = 0. This is quite different from the semi- 
classical phase diagram of Fig. 3, that predicts that the 
AF2 phase should be more stable below C/ < 2.8t for a 
broad range of t' values. This discrepancy is however not 
surprising given that the classical phase diagram was ob- 
tained from the large-[/ expansion that is bound to fail 
in the small-[/ region of the phase diagram. 

One should note that around t' /t = l/\/2 the topology 
of the non-interacting Fermi surface changes, as depicted 
in Fig. 6. This figure suggests why the ordering wave 
vectors take the above mentioned values. 

It is instructive to compare the transition into the Neel 
phase obtained here, with the known Hartree-Fock re- 
sult for the half-filled Hubbard model'^^ (dashed line in 
Fig. 5). The VGA results follow closely the Hartree-Fock 
results for t' < 0.7, whereas for higher levels of frus- 
tration, the non-trivial disordered region is revealed by 
VGA, followed by the (tt, 0) magnetic phase (the latter 
was not considered by the authors of Ref . 37) . 

We note that all the transitions on the phase diagram 
have been found to be first order (with a possible ex- 
ception of very low U values where a reliable solution is 
progressively harder to obtain) . By this we mean that the 
magnetic order parameter jumps as the transition line is 
crossed, so that either no magnetic solution was found in 
the PM phase (typical the low U region) or alternatively, 
since hysteresis is expected, its free energy was found to 
be actually higher than that of the PM solution (typically 
the case for U/t > 6). 

The two AF phases mentioned above, (7r,7r) and (7r,0) 
have been found not only in studies of the J1-J2 Heisen- 
berg model, as already mentioned, but also in work on 
the frustrated Hubbard model by Mizusaki and Imada^^. 
There, the (7r,0) phase is, perhaps confusingly, called a 



"stripe" phase. By using a path-integral Monte-Garlo 
technique, they found a phase diagram very similar to 
ours. In addition, these authors find a third phase with 
a larger periodicity in real space corresponding to the 
ordering vector Q = (7r,7r/2), which exists in a narrow 
region of t'/t between 0.6 and 0.8 for U > 7. This is 
precisely the region where neither (7r,0) nor (7r,7r) phases 
have been found stable in this work. It may well be 
that another magnetically ordered phase (probably in- 
commensurate^^) exists in this intermediate region. Un- 
fortunately, the clusters used in this work, see Fig. 4, 
are not suitable to check for the stability of the (7r,7r/2) 
phase. 

Apart from the possibility of some non-trivial mag- 
netic order, the question remains open whether the 'dis- 
ordered' phase around t/t' — \f2 can actually be formed 
by spin singlets sitting on bonds, as e.g. in the so-called 
valence-bond solid (VBS). This phase is characterised by 
spontaneously broken translational symmetry, but pre- 
serves the spin-rotational symmetry and is considered to 
be the most likely candidate for the intermediate phase 
around the boundary between the two magnetic phases 
of the Ji- J2 Heiscnberg model (see Ref.24 and references 
therein). Another possibility is the spin liquid phase, 
which is similar to VBS, but preserves the translational 
invariance of the system and can be interpreted as a res- 
onating valence bond (RVB) state. 

As it stands, the VGA approach does not permit to 
study the various possibilities for a VBS or spin liquid 
phase. Therefore, at present, we cannot judge whether 
such an order may exist in the paramagnetic region or 
how far it extends into the low-f/ part of the phase dia- 
gram. It is interesting that a quantum spin liquid state 
has been recently shown to exist in the Hubbard model 
on a square lattice at half filling by numerical stud- 
ies on finite clusters^^ using path-integral Monte-Garlo 
mcthod'^^. This state, observed in a narrow region of U jt 
values between 4 and 7, falls between the metallic param- 
agnetic state and the magnetic insulator and is believed 
to be caused by charge fluctuations in the vicinity of the 
Mott transition (see more on that in Sect. VB). Unlike 
in the magnetically ordered state, the quantum spin liq- 
uid is characterized by the absence of any sharp peaks 
in the spin correlation function S'(q). It is certainly an 
intriguing possibility that should be verified with other 
existing numerical methods. 

Naturally, if the VBS or quantum spin liquid solution 
happens to have a lower energy than either of the two 
magnetic phases discussed in this work, this would lead 
to further enlargement of the range of t' values where 
no long-range magnetic order is found. In this sense the 
VGA method only gives a lower bound on the extent 
of the non-magnetic hatched region in Fig. 5, which for 
e.g. Ujt = 9 exists in the range of 0.77 < t'/t < 0.82. 
It is useful to compare these figures with the exact di- 
agonalization results'^^ for the J1-J2 Heisenberg model 
discussed in Sec. II: There, the non-magnetic region ap- 
pears in the range 0.4 < J2/J1 < 0.6. For the Hubbard 
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model, this translates into the large U limit and the win- 
dow 0.63 ^t'/t< 0.78, which is indeed much wider than 
predicted by VGA. 



B. Metal-insulator transition 

We next turn to the subject of the metal-insulator 
transition in the frustrated Hubbard model at half-filling 
in the absence of long-range order. Since there is no bath 
in VGA as we define it, metallic states are less favored 
than in GDMFT. In the normal state, the bath present 
in GDMFT or DMFT can play the role of a metallic or- 
der parameter. While metallic states can occur in VGA, 
they cannot occur as first order transitions because of the 
absence of this metallic order parameter. So, contrary to 
the case of GDMFT^^'''", the Mott transition cannot be 
observed as a cusp or discontinuity in the site double oc- 
cupancy (n^ni) - the dependence of this quantity on U/t 
is a very smooth monotonic curve. 

The Mott transition is however firmly established in 
the half-filled Hubbard model, and can indeed be ob- 
served by analyzing the spectral function A(k, oj). There- 
fore, for the purpose of this study, we define a "metal" as 
a state with non-vanishing spectral function at the Fermi 
level, A(k, iv = 0). We note in passing that the latter def- 
inition is actually broader than saying that there exists 
a well-defined Fermi surface in the ground state, for the 
following reasons. Firstly, the regions of non- vanishing 
A(k, to = 0) need not form a closed surface, but may in- 
stead have a shape of isolated arcs (c.f. the well-known 
Fermi arcs as revealed by ARPES spectroscopy^^ in the 
underdoped cuprates). Secondly, the Landau picture of 
a metal predicts an infinite lifetime for the quasiparti- 
cles at the Fermi surface, equivalent to the requirement 
of delta-function shape for A(k, uj — 0) at the Fermi en- 
ergy. 

Since the presence of a long-range magnetic order 
opens up a gap at the Fermi surface, we intentionally 
suppress the possibility of magnetic ordering. Only the 
effect of short-range magnetic correlations is included in 
the cluster. This is an established practice used to obtain 
the parameters of the metal- insulator transition-'^^'^*'^^''^^. 

The procedure we have adopted is as follows. For a 
given value of t', we plotted the function yl(k, a; = 0) 
across the Brillouin zone (BZ) for several values of the 
interaction U, with the Lorentzian broadening rj = 0.05 
used to account for the imaginary part of A(k, w). The 
point where the spectral function vanishes everywhere 
in the BZ (as U is increased) marks the transition from 
metallic to insulating state. In the present approach, the 
transition appears as a crossover, that strictly speaking 
at zero temperature should be a second-order transition. 
On the anisotropic triangular lattice, it was found that 
as a function of t' , the Mott transition goes from second 
order at small t' to strongly first order at large t' through 
a tricritical point"'^^'"'^. 

Our "crossover" transition line together with error bars 



is plotted in Fig. 5. We see that with increasing frus- 
tration t' , the critical value Uc{t' /t) increases monotoni- 
cally. We note an important difference between the low- 
t' region and that for t' > 0.7t. In the former case of 
almost perfect nesting, the effect of short-range corre- 
lations is strong, leading to a surprisingly low value of 
Uc ~ 2t. At first sight, this is too different from the 
well-known DMFT result^^ for the Mott metal-insulator 
transition, U « 12t. It must be noted however that the 
single-site DMFT approach^'^ does not take into account 
short-range magnetic correlations, as opposed to clus- 
ter methods such as GDMFT or VGA. In addition, to 
study the Mott transition, DMFT assumes the presence 
of large frustration to prevent magnetic long-range or- 
der. Hence, the DMFT result should be compared with 
the region t' > 0.7t in our phase diagram Fig. 5 where 
magnetic order is naturally absent and the critical inter- 
action strength is Uc ~ 5t. In real experiments, that is 
where we believe true Mott physics would be observed. In 
any case, it is physically expected that in two dimensions 
the critical U for the Mott transition depends on frustra- 
tion, as pointed out in GDMFT study of the anisotropic 
triangular lattice^^. 

It is encouraging that values of Uc{t' /t) very similar to 
ours have been found in the path-integral Monte-Garlo 
study^^ of the same model, yielding Uc = 3< at t' /t = 0.25 
and Uc ~ 5t for t'/t = 0.8. However, the authors of 
Ref. 21 have not excluded the possibility of long-range 
magnetic order, which is why they observe the Mott tran- 
sition happening at infinitesimally small values of U in 
the case of perfect nesting t' — 0. We stress that this 
is in complete agreement with our data, although we in- 
terpret this as an opening of the AF gap at the Fermi 
surface rather than Mott transition into a phase with no 
broken spin-rotational symmetry. A somewhat poorer 
agreement is seen with the results of the optimized varia- 
tional Monte Garlo (VMG) method in Ref. 20. There, the 
authors obtain the value Uc ~ 7t in the region \t' /t\ < 0.5 
that they studied. 

Our value of Uc ^ 5t around t'/t = 0.7 should also be 
compared with the recent GDMFT results for strongly 
frustrated lattices, such as the triangular lattice^'^ with 
Uc = 10.5t and the asymmetric square lattice^^ {t' = t 
along only one diagonal) with Uc ~ St. Although 
GDMFT predicts higher U values for the Mott metal- 
insulator transition than VGA in these cases, they clearly 
fall into the same ballpark. However, for the square 
lattice without frustrations {t' = 0) the 4-site cluster 
GDMFT gives a value'^^ for the Mott transition Uc « 5t, 
quite a bit larger than our result Uc ~ 2t. Similar to 
GDMFT, a value Uc ~ 6t was also found in Quantum 
Monte Garlo studies'*^ at t' — 0. This discrepancy is 
however not surprising since it is well known**^, for the 
reasons mentioned at the beginning of this section, that 
VGA method (without a bath) tends to overestimate the 
effect of interactions compared with GDMFT, thereby 
yielding smaller values of Uc for the Mott transition. 
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FIG. 7: (Color online) The phase diagram of the t — t' — U 
Hubbard model obtained with the VGA based on a 8-site 
cluster. The solid lines denote the phase boundaries between 
the (tTjTt) and {n,Q) antiferromagnetic phases and the A^2_y2 
SC phase. The hatched area on the phase diagram (red-and- 
white online) denotes the critical region where neither (tt, tt) 
nor (tt, 0) order could be found. The coexistence region be- 
tween AF and SG phases (shaded area, cyan online) is con- 
tained between two black lines that meet around t' — 0.7t, 
with the (blue) line in between indicating the area where the 
free energies of the would-be separate SG and AF phases be- 
come equal (c.f. the point U = 2A5t in Fig. 8b). Triangles 
and filled circles denote points on the phase boundary where 
an order parameters sustains a discontinuity at a first-order 
phase transition; short dashes inside the coexistence region 
mark the points where total energies of the AF and SG phases 
are equal. 



VI. INTERPLAY BETWEEN MAGNETISM 
AND SUPERCONDUCTIVITY 

A. Results and discussion 

Our final VGA phase diagram is shown in Fig. 7 for the 
largest 8-site cluster studied. Most interestingly, in addi- 
tion to the two magnetic phases discussed in the previous 
section, a d-wave superconducting solution (with dx^~y^ 
symmetry) comes out naturally from the VGA calcula- 
tions for low values of U /t in the phase diagram. It is 
clear from Fig. 7 that the frustration tends to destroy the 
Neel phase and stabilise the SC solution. In the whole 
range \t'/t\ < 1, we did not find any superconducting 
regions with stable dxy symmetry of the order parame- 
ter, although there are indications^^ that such a phase 
would become stable in the case of (unrealistically large) 
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FIG. 8: (Color online) Details of the AF and SC phase coex- 
istence for a 2x3 cluster at t' = 0.2t. (a) Expectation values 
of SG {D, circles and red line) and AF {M, squares and black 
line) order parameters as a function of increasing interaction 
U; note different scales for the two quantities plotted. Values 
Ui and U2 denote the positions of the first-order phase tran- 
sitions into/from the coexistence phase, (b) Blow up of the 
free energy F — fl + fJ,{n) as a function of U near U = 2.5t is 
shown for all three phases studied: AF phase (blue solid line), 
SG phase (red dash-dotted line) and for the coexistence phase 
where both D and M order parameters are non-zero (dashed 
black line). The latter phase has lower energy than the other 
two in the whole coexistence region of 2.3 < U/t < 3.9. 

frustration strength t' /t > 1.1. 

In the low-C/ region of the phase diagram, both AF and 
d-wave SG solutions are stable, hence the one with lowest 
free energy F would win. Since we work at a fixed particle 
density n (half- filling) , we perform the Legendre trans- 
form to obtain the free energy from the grand-canonical 
potential f2: 

F = E -TS + fi{n). (17) 
Moreover, since the calculations are done at zero temper- 
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FIG. 9: (Color online) The expectation values of SC {D, cir- 
cles and red line, left-hand scale) and two types of AF order 
parameters (shown with squares and black lines), calculated 
on a 2x3 cluster for (a) t' — 0.5t, the AF (tt, tt) phase order 
parameter Af (right-hand scale), and (b) t' = 0.8t, the (7r,0) 
phase order parameter Q (right-hand scale). Unlike in (a), 
no coexistence region is seen in case (b) , where the first-order 
transition occurs at Uc = 4.7f. 



ature. the Gibbs free energy in this case is identical to the 
total energy E of the system. We find that the transition 
between the magnetic and the SC phases is first order, 
with the blue line in Fig. 7 (in the centre of the shaded re- 
gion) denoting the points where the total energies of the 
two competing phases (no coexistence allowed) become 
equal. The free energy of both phases is illustrated in 
more detail by the dot-dashed and solid lines in Fig. 8b. 

Most interestingly, the lowest energy solution (dashed 
line in Fig. 8b) corresponds to a new phase, shown as 
shaded area in Fig. 7, where both magnetic and SC or- 
der parameters are non-zero. This is a phase with a true 
homogeneous coexistence of the magnetic and supercon- 
ducting phases, which one may want to call an antifer- 
romagnetic superconductor to emphasize the difference 
from a more usual inhomogeneous coexistence observed 
e.g., at a first-order transition. 

The details of the transition between the coexistence 
phase and the pure AF and SC phases are illustrated in 
Figs. 8a and Fig. 9a that show the dependence of the 
corresponding order parameters, A'l and D, on the inter- 
action strength U. AsU decreases below Ui in Fig. 8a, we 
first observe a first-order transition from the coexistence 
phase into a pure d-wave SC state, where the AF order 
parameter plunges to zero and the SC order parameter 
sustains an upward jump as the coexistence phase ceases 
to exist. As the interaction increases above U = U2, 
there is a similar transition from the coexistence phase 
into the pure antiferromagnct (7r,7r), although this time 
the transition appears more continuous (sec Fig. 8a). 

Clearly, our phase diagram Fig. 7 shows that frustra- 
tion t' /t favors the SC phase as long as it is not too large. 
Indeed SC becomes more stable and occupies a broader 
region of the phase diagram as frustration increases at 
low to intermediate interaction strength U, until t' /t be- 
comes large enough for the AF2 (tt, 0) phase to decrease 
the area occupied by the SC phase. The latter transition 
is of first order and is accompanied by a sharp jump in 
the values of the respective order parameters a,t U = Uc, 
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FIG. 10: (Color online) The phase diagram of the t — t' — U 
Hubbard model obtained with the VGA based on an (a) 4- 
site cluster and (b) 6-site cluster. The hatched area (red-and- 
white online) denotes the critical region where neither (tt, tt) 
nor (tt, O) order could be found or when the paramagnetic 
solution had lower energy then any of the AF phases. The 
coexistence region between SC anf AF phases (shaded area, 
cyan online) is shown. The solid lines, as well as open trian- 
gles and short dashes on the phase boundaries have the same 
meaning as in Fig. 7. 



as one can verify on Fig. 9b for t'/t = 0.8 for the 6 
site cluster. Sometimes a very narrow hysteresis region 
(AU/t < 0.2) has been observed, depending on whether 
the transition is approached from above or below the crit- 
ical value of Uc- 

In the region of t'/t « 0.75 of the phase diagram Fig. 7, 
the SC phase has a direct boundary with the disordered 
non-magnetic phase as a function of increasing U. Given 
the possibility for the existence of the VBS order in that 
region, this opens up the interesting possibility of a di- 
rect transition from VBS into the SC state. Since the 
two phases have symmetry groups that are not related 
to each other, according to Landau theory they would be 
separated either by a first order transition or by a coex- 
istence phase. Or, beyond the Landau paradigm, they 
could lead to an example of deconfined quantum critical- 
ity*^. A similar deconfined critical transition could also 
possibly occur between the VBS and antifcrromagnetic 
phases. 

To assess convergence towards the thermodynamic 
limit, it is useful to compare the phase diagrams ob- 
tained for different cluster sizes. Figure 10 shows the 
VGA phase diagrams obtained from the 4-site and 6-site 
clusters. We see that the main conclusions, namely that 
the d-wave SC phase is the ground state for low U values, 
remains unchanged. However the position of the AFl-SC 
phase boundary shifts to lower U values with increasing 
cluster size. This is to be expected since perfect nesting 
at = and half-filling should lead to antiferromag- 
netism as the true ground state of the Hubbard model at 
arbitrary small values of U. This suggests that perhaps 
the coexistence of SC and AFl at small t' ^ should 
disappear altogether in the thermodynamic limit of in- 
finitely large cluster size. We do see from Figs 10 and 7 
a decrease in the size of the coexistence region as the 
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cluster size increases. 

We also note that the phase boundary of the coUinear 
(tt, 0) phase occurs at lower values of U with increasing 
cluster sizes. For example, comparing Figs. 10 and 7, 
we see that at t'/t = 0.96 the transition from the AF2 
phase into the SC phase occurs at U/t « 3.5 for the 8-site 
cluster instead oiU /t~ 5 that one observes for the small- 
est cluster studied. As in the case of the largest cluster 
studied (phase diagram in Fig. 7 above), all transitions 
between the different phases, including the coexistence 
phase between AF and SC orders, were found to be first 
order. 

We now compare our phase diagram in Fig. 7 with 
that found by other methods. Recently, the path-integral 
Monte Carlo study of the same model by Mizusaki and 
Imada^^ has revealed a phase diagram where the mag- 
netic and PM regions are in very good agreement with 
our Fig. 5, apart from an extra magnetic phase that the 
authors of Ref. 21 observe between the AF and the (tt, 0) 
phases. This however does not contradict our results 
since even more complicated phases, such as incommen- 
surate magnetic order, may be possible but are beyond 
reach of quantum cluster methods such as VCA. Another 
important difference is the existence of the quantum spin 
liquid phase found in Ref. 21 that we commented upon 
in Sect. VA. The possibility of superconducting phase 
has not been addressed however in Ref. 21. This would 
have been very instructive in light of our findings. 

Another recent study of the half-filled Hubbard model 
has been performed recently by Yokoyama et alP'^ using 
optimization VMC. They considered both the (tt, tt) AF 
phase and dx'^-y2 superconductivity. It is puzzling that 
the AF phase has been found to be limited to the region 
\t' /t\ < 0.2 only, in contradiction with both our work 
and the known results of the J1-J2 model in the strong 
coupling limit. As the authors themselves suggest, this 
discrepancy is most likely due to a bad choice of the vari- 
ational state and/or the limitations of the VMC method. 
In that work, the d-wave SC state is found to be most 
stable in the vicinity of the Mott transition {Uc = Q.7t) 
for a narrow range of 0.2 < \t'/t\ < 0.4, although the 
authors mention that a small magnitude of the SC gap 
survives often to very small values of U/t, which would 
be in agreement with our results showing that SC exists 
in the whole region of < \t' /t\ < 1 down to C/ = 0. 
Unfortunately, the limited range of \t'/t\ < 0.5 studied 
in Ref. 20 does not help to shed light on the existence of 
other AF orderings, such as the (tt, 0) state, or indeed on 
the long sought after quantum spin liquid state, claimed 
to have been observed convincingly in Ref. 21. 



B. Comparison with results on the anisotropic 
triangular lattice 

To conclude this section, we contrast the results with 
those obtained on the anisotropic triangular lattice^*' 
First of all, on that lattice CDMFT shows that phase 



transitions between ordered phases occur at the same 
location as the Mott transition when the latter is first 
order^^. This does not happen when the transition is 
second order^^. In our case, the Mott transition is al- 
ways second order and the transitions between ordered 
phases do not coincide with the Mott line as we can see 
by comparing the line with error bars on Fig. 5 with the 
phase diagram in Fig. 7. 

On the anisotropic triangular lattice, the transition be- 
tween d-wave superconductivity and antiferromagnetism 
is always first order with no coexistence region, in con- 
trast with our case where (tt, tt) antiferromagnetism is 
separated from d-wave superconductivity by a coexis- 
tence phase. At larger frustration however, d-wave super- 
conductivity is separated from (tt, 0) antiferromagnetism 
by a first order transition. Since the triangular lattice 
has geometric frustration, not only frustration induced 
by interactions, one can speculate that it is the larger 
frustration on the anisotropic triangular that leads to 
the disappearance of the coexistence phase. However, 
a.t t' — both problems become identical and there is 
a clear disagreement between the results that can come 
only from differences in the two calculational approaches, 
CDMFT vs VCA. 

It is rather striking that on the anisotropic triangu- 
lar lattice, the d-wave order parameter is largest as a 
function of U jt when it touches the first order boundary 
with the antiferromagnetic phase. As can be seen from 
Figs. 8a and 9a, this also occurs in our case when d-wave 
superconductivity touches the homogeneous-coexistence 
phase boundary. Also, the maximum value of the order 
parameter increases in going from t' = 0.2t to t' = 0.5t, 
as in the case of the anisotropic triangular lattice. At 
larger frustration, the trend as a function of t' jt reverses 
on the latter lattice. In our case, we see that at large 
frustration t' = on Fig. 9b, the d-wave order param- 
eter reaches its maximum value before it hits the first 
order boundary with the (tt, 0) phase. 

It would clearly be interesting to compare what are 
the predictions for our case of other quantum cluster 
approaches, such as CDMFT and DCA. The weak to 
intermediate coupling Two-Particle Self-Consistent Ap- 
proach'^^ suggests that at small values of U superconduc- 
tivity disappears in favor of a metallic phase and that 
at t' = antiferromagnetism dominates. We expect the 
predictions of VCA to be more reliable at strong cou- 
pling. 



VII. SUMMARY AND CONCLUSIONS 

We have analyzed the phase diagram of the half-filled 
Hubbard model as a function of frustration t'/t and in- 
teraction U /t using both analytical and numerical tech- 
niques. 

The classical analysis based on the resulting large-t/ ef- 
fective spin Hamiltonian allowed us to draw the classical 
magnetic phase diagram shown in Fig. 3. Of course, this 
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classical approach remains completely oblivious to the 
role of quantum fluctuations and in addition is designed 
for the high-C/ segment of the phase diagram where elec- 
trons are localized. 

To treat quantum fluctuations as well as the possibil- 
ity of delocalization, we used the Variational Cluster Ap- 
proximation. Because of the nature of the VGA method, 
it relics on the choice of a finite (necessarily small) clus- 
ter on which the problem can be solved exactly. In order 
to the finite-size effects, we have analyzed the results for 
clusters of 4, 6 and 8 sites (the latter being almost at 
the limit of what can be achieved with today's power- 
ful supercomputers using the exact diagonalization algo- 
rithm). Although the exponentially increasing computa- 
tional cost did not permit us to study larger clusters, the 
apparent similarities between the 6- and 8-sitc cluster so- 
lutions allow us to conclude that the VGA calculations 
reported in this study arc close to convergence with re- 
spect to increasing cluster size. 

Independent of the cluster size used, the key features 
of the resulting phase diagram are as follows. At large 
values of U the VGA results agree qualitatively with 
the classical \aige-U expansion (Sec. Ill) of the Hubbard 
Hamiltonian and with the known results for the related 
J1-J2 Heisenberg spin model (Sec. II), which show that 
two competing AF phases AFl and AF2 with ordering 
wave-vectors Qi = (7r,7r) and Q2 = (7r,0) exist for frus- 
trations lower and higher, respectively, than some critical 
value given roug hly by t'Jt = l/^/2 « 0.71. This is where 
the Fermi surface changes topology in the non-interacting 
case (see Fig. 6). These two magnetic phases are sepa- 
rated by a disordered region where we don't find non-zero 
values of either order parameter. The Heisenberg model 
studies point to possible existence of an exotic valence 
bond sohd (VBS) phase around the critical frustration 
value t'^. Although the direct study of the VBS phase re- 
mains beyond reach of quantum cluster methods such as 
VGA, it is encouraging that the obtained phase diagram 
exhibits a region around the critical frustration tc where 
no long-range AF order could be found. 

We have also addressed the issue of the role of frustra- 
tion on the metal-insulator transition that is known to 
exist in the Hubbard model at half-filling. The distinc- 
tion between metallic and insulating phases was based 
on the analysis of the spectral function yl(k, w) at the 
Fermi level, which is vanishing in the insulating ground 
state. We find that the value of the interaction strength 
Uc where the insulator appears rises monotonically as a 
function of frustration strength t'. The value of Uc turns 
out to be surprisingly low for small t' values {Uc ~ 2t). 
Low values are expected from the effect of short-range an- 
tiferromagnetic fluctuations that are particularly strong 
near perfect nesting at t' — 0. Nevertheless, comparisons 
with other results definitely suggest that VGA overem- 
phasizes the effect of U so that the insulator-metal tran- 
sition at half-filling should be closer to the GDMFT and 
QMC values Uc ^ 5 — 6t. For large frustration near 
t'^/t = ~ 0.71, the transition occurs at larger val- 



ues of Uc « 5t. Since neither the AFl or AF2 phases 
are stable in this highly frustrated region, this is where 
one would experimentally be more likely to see a gen- 
uine Mott transition at finite temperature where ordered 
phases are absent. 

Most importantly, the VGA method allowed us to 
study the region of the phase diagram with low to inter- 
mediate interaction strength, which is inaccessible in the 
large-i7 expansion or the J1-J2 Heisenberg spin model. 
We find that even at half-filling, where the tendency 
towards the antiferromagnetic ordering is strong, frus- 
tration allows d-wave superconductivity to appear for a 
range of values of U /t that generally increases with frus- 
tration since the latter is detrimental to (tt, tt) antiferro- 
magnetism. With frustration in the range t'c/t « 1/v^, 
both (tt, tt) and (tt, 0) antiferromagnetism disappear but 
d-wave superconductivity survives for Uc ^ 5t. Increas- 
ing frustration further favors the (tt, 0) phase but d-wave 
superconductivity continues to appear at smaller values 
of U/t. All the phase transitions are first order, except 
possibly the transition from the coexistence phase into 
the (tt, tt) magnet, which is weakly first order (note a 
very small jump in the value of SC order parameter at 
the U = U2 phase boundary in Fig. 8 a). In the case of 
(tt, tt) antiferromagnetism, the transition to pure d-wave 
superconductivity occurs through a region where both 
phases coexist homogeneously. Finite size analysis sug- 
gests that this coexistence region is relatively robust, al- 
though its boundaries shrink with increasing cluster size. 
Goexistence may disappear in the thermodynamic limit. 

An important prediction of our study for experiments 
is that d-wave superconductivity may appear by applying 
sufficiently high pressure on the half-filled parent com- 
pounds of high-temperature superconductors. This type 
of transition is observed in layered BEDT organics^ and 
can be explained by the Hubbard model^^'^^. Hence pos- 
itive results of such an experiment on the cupratcs would 
spectacularly help to establish definitively the electronic 
origin of d-wave superconductivity. 

It would be interesting to pursue the issues addressed 
in this work with other quantum cluster approaches and 
also to study the case of doped Hubbard model away from 
half-filling, with possible comparison with the results of 
the much studied t-J model. Anticipating on the results, 
it should be easier to reach the d-wave superconducting 
state by applying pressure on a slightly doped insulating 
parent than on the half-filled insulator. These issues, of 
much relevance to the physics of high-tcmpcrature super- 
conductivity in the cuprates, are left for future studies. 
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